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Abstract 

Research  results  for  AFOSR  grant  F49620-03- 1  -02 1 2  for  the  period  1  April,  2003  to  31 
December,  2006  is  described.  The  objectives  of  this  research  were  to  develop  and 
validate  classical  force  fields  for  different  model  compounds  representative  of  a  range  of 
ionic  liquids;  compute  a  wide  range  of  physical  properties  for  these  model  compounds; 
obtain  molecular-insight  into  how  property  variations  are  related  to  structure;  and 
develop  new  simulation  tools  to  accurately  compute  melting  points  and  gas  solubilities  in 
ionic  liquids.  All  of  these  objectives  were  achieved.  Force  fields  for  a  range  of 
imidazolium-  pyridinium-  and  triazolium-based  ionic  liquids  were  developed  and 
published.  Properties  including  densities,  heat  capacities,  cohesive  energy  densities, 
enthalpies  of  vaporization,  gas  solubilities  and  diffusion  coefficients  were  computed. 
Liquid  structure  and  these  properties  were  relates  to  molecular  interactions  identified  in 
the  simulations.  A  new  rigorous  melting  point  prediction  simulation  method  was 
developed  and  applied  to  test  systems  including  NaCl,  benzene  and  triazoic.  It  is 
currently  being  applied  to  ionic  liquid  systems.  A  new  semi-grand  ensemble  simulation 
method  for  predicting  liquid-liquid  equilibrium  and  vapor-liquid  equilibrium  was 
developed  and  published. 
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1.  Detailed  Objectives 


Three  major  objectives  were  stated  in  the  original  proposal:  1)  Develop  and  validate 
molecular  mechanics  force  fields  for  different  model  compounds  representative  of  a 
range  of  energetic  ionic  liquids;  2)  Compute  a  wide  range  of  physical  properties  for  the 
model  compounds  and  related  species,  including  densities,  isothermal  compressibilities, 
coefficients  of  thermal  expansion,  surface  tensions,  viscosities,  diffusivities,  melting 
points  and  heat  capacities;  3)  Obtain  molecular-level  insight  into  how  property  variations 
are  related  to  structure,  and  share  this  information  with  synthetic  groups  to  help  guide 
their  efforts.  A  fourth  objective  was  added  soon  after  the  start  of  the  project:  4)  Develop 
new  simulation  tools  to  accurately  compute  melting  points  and  gas  solubilities  in  ionic 
liquids. 

2.  Status  of  Effort 

Objective  1:  New  all-atom  forcefields  for  imidazolium-based,  alkylpyridinium-based  and 
triazolium- based  ionic  liquids  were  developed  and  published  in  the  open  literature. 
Liquid  and  solid  properties  of  these  compounds  were  computed  and  compared  with 
experimental  data  provided  by  collaborators  to  obtain  experimental  validation. 

Objective  2:  A  range  of  thermodynamic  and  transport  properties  for  pure  pyridinium-, 
triazolium-,  and  imidazolium-based  ionic  liquids  were  computed. 

Objective  3:  Based  on  the  results  above,  the  mechanism  for  high  CO2  solubility  in  ionic 
liquids  was  provided.  In  addition,  insight  into  melting  point  trends  was  provided  by 
showing  that  high  cohesive  energy  densities  correlate  with  high  melting  point.  Unusual 
experimental  observations  made  regarding  liquid-liquid  phase  behavior  between  ionic 
liquids  and  alcohols  was  rationalized  in  terms  of  strength  of  interactions  between  the 
cation,  anion  and  alcohol. 

Objective  4:  A  new  free  energy-based  simulation  method  for  rigorously  computing 
melting  points  of  atomic  species  was  developed  and  extended  to  molecular  systems.  Two 
new  semi-grand  ensemble  simulation  methods  were  developed  specifically  for  computing 
gas  solubility  in  liquids.  This  approach  is  currently  being  applied  to  directly  compute 
isotherms  for  ionic  liquids. 
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3.  Accomplishments  /  New  Findings 


3.1  Development  of  a  forcefield  for  heterocyclic-based  energetic  ionic 
liquids  formed  from  triazolium  precursors 

C.  Cadena  and  E.  J.  Maginn,  "Molecular  simulation  study  of  some  thermophysical  and 
transport  properties  of  triazolium-based  ionic  liquids",  Journal  of  Physical  Chemistry  B , 
2006,  110,  18026-18039 

David  M.  Eike  and  Edward  J.  Maginn,  “Atomistic  Simulation  of  Solid-Liquid 
Coexistence  for  Molecular  Systems:  Application  to  Triazole  and  Benzene”,  Journal  of 
Chemical  Physics,  2006,  124,  164503. 

3.1.1  Motivation 

Drake  and  co-workers 1  recently  announced  the  synthesis  of  three  new  families  of 
heterocyclic-based  energetic  salts.  The  cations  for  these  salts  are  1-H- 1,2,4  triazolium, 
4-amino- 3,2,4  triazolium  and  1-H- 1,2,3  triazolium.  These  cations  arc  shown  in  Fig.  1. 


Figure  l:  The  l-H-1,2,4  triazolium,  4-amino- 1,2,4  triazolium  and  I -H- 1,2,3  triazolium 
salts  synthesized  by  Drake  and  co-workers.  X  can  be  NO},  CIO 4 ,  or  N(NOi)i. 

These  salts  have  a  number  of  favorable  properties,  including  good  thermal  stability, 
relatively  low  melting  points,  and  high  energy  content.  Although  the  crystalline  densities 
and  some  melting  points  have  been  measured  in  the  lab,  liquid  phase  properties  have  not 
been  measured  due  to  the  difficulty  of  the  experiments2.  For  this  reason,  we  developed  a 
forcefield  for  the  nitrate  and  perchlorate  forms  of  these  salts  and  compute  liquid  phase 
properties.  We  also  developed  a  forcefield  for  the  l-methyl-4-amino-l,2,4-triazolium 
cation,  as  this  compound  has  recently  been  made  in  Prof.  Shreeve’s  group. 

3.1.2  Functional  form  of  the  forcefield 


The  functional  form  of  the  forcefield  used  in  this  work  was  similar  to  that  used  in  our 
earlier  work  with  iinidazolium-based  ionic  liquids.  The  total  energy  of  the  system  Vwi  is 
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where  the  individual  terms  account  for  bond  stretching,  bond  bending,  torsional  rotation, 
improper  angle  bending,  repulsion-dispersion  interactions  and  electrostatic  interactions. 

The  parameters  for  the  forcefleld  were  determined  by  first  focusing  on  the  molecular 
liquid  1-H- 1,2, 4-triazole.  This  was  done  because  it  is  a  representative  compound  for  the 
all  the  triazolium  cations  and  there  are  experimental  crystal  structure  and  liquid  density 
data  available.  In  addition  to  the  triazolium-based  materials,  a  forcefield  for  l-n-butyl-3- 
methylimidazolium  nitrate  ([bmim][NOj])  was  also  developed  since  we  have  experience 
with  this  cation  as  well  as  experimental  data  for  the  ionic  liquid  itself.  This  allowed  us  to 
assess  the  accuracy  of  the  NOj  parameters. 

To  determine  optimized  geometric  structures,  ab  initio  calculations  were  run  at  the 
B3LYP/6-3 1 1+G*  level  of  theory.  Partial  charges  of  all  of  the  structures  were  determined 
by  applying  the  CHELPG  algorithm.  Cations  and  anions  were  simulated  separately  and 
the  net  charge  was  forced  to  be  +1  /  -1.  Bond  and  bond  angle  force  constants  outside  of 
the  triazole  ring  were  computed  by  perturbing  the  corresponding  minimum  energy  value 
and  then  fitting  the  energy  difference  to  the  appropriate  harmonic  expression.  The  same 
procedure  was  followed  for  computing  the  bond  angle  force  constants  inside  of  the  ring 
(5  bonds  and  5  angles).  Dihedral  parameters  involving  atoms  in  the  ring  were  set  to  zero. 
Improper  torsion  parameters  were  computed  from  the  frequency  modes  analysis  through 
the  visual  identification  of  the  vibrational  modes.  Lennard-Jones  parameters  for  the  heavy 
atoms  were  estimated  using  values  for  similar  atoms  in  the  CHARMM  forcefield 
database,  For  the  hydrogens,  the  most  effective  set  of  Lennard-Jones  parameters  was 
determined  by  computing  densities  for  a  range  of  different  parameter  values  and 
comparing  predicted  densities  against  experiment. 


To  assess  the  accuracy  of  the  intramolecular  parameters,  calculations  were  first 
performed  on  l-H-l,2,4-triazole.  Figure  1  shows  the  minimized  structure  of  this 

molecule.  The  crystal  structure  of  this  compound  is  known, 
and  atomic  coordinates  are  available  through  X-ray 
diffraction3,  low-temperature  electron  diffraction4, 
microwave  data5,  and  electron  diffraction  studies6.  Table  1 
shows  a  comparison  between  the  ab  initio  calculations  (AI) 
and  the  experimental  data.  As  can  be  seen,  agreement  is 
excellent  and  is  certainly  within  the  experimental 
uncertainty  and  fluctuations. 


Figure  1:  Structure  of  1- 
H- 1 ,2,,4-triazole 


The  cation  4-amino- 1, 2, 4-triazolium  required  additional 
dihedral  torsion  and  Lennard  Jones  parameters,  which  were 
taken  from  similar  compounds  in  the  CHARMM  database. 
All  of  the  force  constants  and  the  Lennard-Jones  parameters  for  the  [bmim]  cation  were 
taken  from  our  previous  work7.  Partial  charges  were  recomputed  through  ab  initio 
calculations  so  that  a  formal  charge  equal  to  unity  was  obtained.  A  complete  listing  of  the 
parameters  and  atom  definitions  was  given  in  'fables  1-7  in  the  Supplemental  Material 
section  of  the  2004  report. 
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Table  1:  Computed  and  experimental  bond  lengths  and  angles  for  J-H- 1.2, 4-triazole. 
Al=  ab  initio  calculations;  RT,  LT  =  room  temperature  and  low  temperature  X-Ray 
diffractions,  respectively ;  MW  -  microwave;  ED=electron  diffractions. 


Bonds,  A 

AI 

KT3 

LT4 

MW5 

ED6 

Exp 

average 

%  dev. 

N2-C3 

1.321 

1.330 

1.323 

1.328 

1.329 

1.328 

-0.53 

N2-N1 

1.355 

1.354 

1.359 

1.381 

1.380 

1.369 

-1.02 

C3-N4 

1.363 

1.353 

1.359 

1.354 

1.348 

1.354 

0.66 

C3-H7 

1.080 

- 

0.930 

1.078 

1.054 

1.021 

5.78 

N4-C5 

1.318 

1.352 

1.324 

1.280 

1.305 

1.319 

-0.08 

C5-N1 

1.350 

1.344 

1.331 

1.375 

1.377 

1.357 

-0.51 

C5-H8 

1.080 

- 

0.930 

1.078 

1.054 

1.021 

5.78 

N1-H6 

1.008 

- 

1.030 

0.998 

0.990 

1.006 

0.20 

Angles,  deg 

C3-N2-N1 

101.9 

101.8 

102.1 

102.7 

102.7 

102.3 

-0.39 

N2-C3-N4 

115.2 

114.5 

114.6 

113.0 

113.8 

114.0 

1.05 

N2-C3-H7 

121.6 

- 

114.0 

128.5 

127.0 

121.3 

0.25 

N4-C3-H7 

123.2 

- 

132.0 

118.5 

119.2 

123.2 

0.00 

C3-N4-C5 

102.8 

104.3 

103.0 

106.8 

105.7 

104.7 

-1.81 

N4-C5-N1 

109.8 

107.1 

110.1 

109.0 

108.7 

108.7 

1.01 

N4-C5-H8 

126.6 

- 

127.0 

120.5 

120.3 

123.8 

2.26 

N1-C5-H8 

123.6 

- 

123.0 

130.5 

131.0 

128.2 

-3.59 

N2-N1-C5 

110.3 

112.2 

110.2 

108.5 

108.9 

110.0 

0.27 

N2-N1-H6 

119.9 

- 

124.0 

127.5 

110.9 

120.8 

0.74 

C5-N1-H6 

129.8 

- 

126.0 

124.0 

140.2 

3.1.3  Summary  of  Major  Findings 

Molecular  dynamics  simulations  were  carried  out  with  this  forceficld  using  procedures 
we  have  described  in  our  previous  papers.  The  following  cations  were  simulated:  1,2,4- 
triazolium,  1 ,2,3-triazolium,  4-amino- 1, 2, 4-triazolium,  l-mcthyl-4amino-l,2,4-triazolium 
and  l-rt-butyl-3-methylimidazolium  nitrate.  Anions  include  nitrate  and  perchlorate.  The 
minimized  structures  of  these  compounds  determined  from  the  DFT  calculations  are 
depicted  in  Figure  2.  1-H-I,2,4  triazole  was  also  simulated  in  the  crystalline  and  liquid 
state. 

Table  2  lists  the  computed  and  experimental  densities  for  l-H-l,2,4-triazo!e.  As  can  be 
seen,  the  overall  densities  agree  reasonably  well  with  the  experiments,  but  the  density  of 
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the  crystalline  phase  in  particular  is  too  low.  We  believe  this  is  due  to  inaccurate 
treatment  of  the  hydrogen  bonding  that  takes  place  in  the  crystal.  This  also  leads  to  an 
under- prediction  of  the  melting  point  for  this  compound  (not  shown).  The  liquid  phase 
densities  are  very  close  to  the  experimental  values  (around  1%),  suggesting  that  the 
forcefield  does  a  better  job  in  the  liquid  phase. 


Figure  2a)  1,2,4-triazoIium  [ I24t]  Figure  2b)  1,2,3-triazolium  [I23tj 


Figure  2c)  4-amino- 1, 2,4- 
triazolium  [4amt] 


Figure  2d)  I -methyl-4-amino- 

1 , 2, 4-triazolium  [ Im4amt]  Figure  2e)  nitrate  [NO 3] 


Figure  2J)  perchlorate  [CIO,} 


Table  2:  Experimental  and  computed 
densities  for  1  -H- 1, 2, 4-triazole 


T  (K)  p*ipt  (g/cm3) 

Pstm  (g/cm3) 

113 

1.448 

1.421 

298 

1.394 

1.346 

426 

1.132 

1.146 

478 

1.089 

1.101 

The  liquid  structure  of  the  triazolium-based  ionic  liquids  was  studied  through  radial 
distribution  functions  (RDFs).  As  an  example,  Figures  3a  and  3b  show  the  center-of-mass 
RDFs  for  the  perchlorate  anions  with  1,2,4-triazoIium  and  4-amino- 1, 2, 4-triazolium, 
respectively.  The  First  peak  is  much  sharper  for  the  1,2,4-triazoIium  cation  than  for  the  4- 
amino- 1,2, 4-triazolium  case.  Part  of  the  difference  can  be  attributed  to  the  larger  cation 
making  the  center-of-mass  more  diffuse.  However,  an  additional  factor  seems  to  be  that 
the  amino  group  weakens  the  interaction  between  the  cation  and  the  anion,  which  is 
consistent  with  its  greater  molar  volume. 


7 


Figure  3a:  COM  RDF  for  1,2,4-triazolium  and  perchlorate.  Solid  line:  +/-,  dashed  line: 
+/+,  dot-dash: 


Figure  3b:  COM  RDF  for  4-amino- 1,2,4-triazolium  perchlorate.  Solid  line:  +/-,  dashed 
line:  +/+,  dot-dash: 

To  better  understand  how  the  anion  organizes  about  the  cation,  we  computed  RDFs  for 
the  cation  center-of-mass  and  the  oxygen  atoms  of  the  anion  as  well  as  the  nitrogen  atom 
(for  NO3  )  or  the  chlorine  atom  (for  CIO4*).  Figure  4  shows  the  result  of  these  calculations 
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for  all  the  ionic  liquids.  Solid  lines  correspond  to  cation  -  O  RDFs,  while  dashed  lines 
are  for  cation  -  N/CI. 


Figure  4:  O-cation  (solid  lines)  and  N/Cl-cation  RDFs  for  all  the  compounds  examined. 


The  results  in  Figure  4  along  with  additional  geometric  analyses  indicate  that  the  plane  of 
the  NO3'  anion  is  oriented  perpendicular  to  the  cation  ring  with  two  oxygen  atoms  facing 
the  ring.  This  makes  sense  from  an  energetic  standpoint,  since  the  oxygen  atoms  have  the 
greatest  negative  charge.  The  tetrahedral  perchlorate  anion  aligns  such  that  three  of  its 
oxygen  atoms  face  the  cation. 
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Densities  for  the  liquid  phase  of  all  compounds  were  computed  as  a  function  of 
temperature.  Unfortunately,  we  do  not  have  any  experimental  liquid  densities  against 
which  to  compare,  so  these  calculations  are  all  predictions.  Table  3  lists  these  results, 
along  with  experimental  melting  points  and  crystalline  densities. 


Table  3:  Computed  liquid  densities  for  various  triazolium-based  ionic  liquids.  No 
experimental  data  is  available,  so  crystalline  density  and  melting  point  is  shown. 


compound 

Tm 

ps(gW) 

373  K 

398  K 

423  K 

448  K 

473  K 

[123t][NO,] 

383 

1.57 

1.609 

1.595 

1.579 

1.561 

1.545 

[1 24t][NOj] 

410 

1.64 

1.634 

1.601 

1.585 

1.568 

1.551 

[4amt][N03j 

342 

1.60 

1.587 

1.577 

1.561 

1.546 

1.529 

ri23t][C10,l 

346 

1.79 

1.738 

1.724 

1.709 

1.693 

1.677 

[124t][C104] 

362 

1.96 

1.748 

1.733 

1.719 

1.706 

1.686 

[4amt][CI04] 

357 

1.81 

1.710 

1.694 

1.681 

1.661 

1.645 

By  examining  the  temperature  and  pressure  dependence  of  the  densities,  volumetric 
expansivities  and  isothermal  compressibilites  may  be  computed.  Heat  capacities  can  also 
be  obtained  from  the  change  in  internal  energy  and  enthalpy  with  temperature.  The 
computed  values  for  these  properties  are  shown  in  Table  4. 


Table  4:  Computed  volumetric  expansivities,  isothermal  compressibilities  and  heat 
capacities  for  six  energetic  ionic  liquids. 


compound 

ap  (K'‘)x  104 

Kt  (bar'1)  xI0s 

Cp  (J/mol  K) 

Cv  (J/mol  K) 

[123tl[N03l 

4.13 

1.17 

398.1 

346.7 

[124tlfN03l 

4.01 

1.12 

403.9 

353.1 

[4amt][N03] 

3.78 

1.08 

463.6 

410.8 

ri23t][C104] 

3.58 

1.23 

409.1 

365.0 

ri24tl[CIO<l 

3.50 

1.14 

409.1 

364.1 

[4amt][CI04] 

3.89 

1.18 

479.6 

419.9 

The  data  in  Table  4  show  that  compounds  are  all  expected  to  have  similar  volumetric 
changes  with  temperature  and  pressure.  The  biggest  difference  is  observed  in  the  heat 
capacities.  The  presence  of  the  amino  group  at  the  4  position  on  the  triazolium  ring 
significantly  increases  the  heat  capacity.  The  other  interesting  observation  is  that  all  of 
the  heat  capacities  are  much  smaller  than  those  seen  with  the  alkyl-imidazolium  species. 
For  example,  the  heat  capacity  of  l-rc-butyl-3-methylimidazolium  nitrate  is  about  a  factor 
of  two  larger  than  any  of  the  triazolium  heat  capacities.  Part  of  this  is  likely  due  to  the 
long  alkyl  group  which  provides  an  additional  mechanism  for  energy  storage. 
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3.1.4  Relationship  to  the  Goals  of  the  Project 

Developing  and  validating  a  transferable  foreefield  for  triazolium  (and  ultimately, 
tetrazolium)  ionic  liquids  is  essential  in  our  ability  to  predict  liquid  properties.  Such 
properties  are  needed  for  engineering  calculations.  Experiments  to  measure  many  of  these 
properties  are  difficult,  and  in  fact  hardly  any  liquid  phase  properties  are  known.  The 
simulation  results  we  have  reported  are  the  very  first  examination  of  this  class  of  ionic 
liquids,  and  provide  much  needed  physical  property  estimates. 

3.2  Development  and  validation  of  a  force  field  for  pyridinium-based 
ionic  liquids 

Cesar  Cadena,  Qi.  Zhao,  Randall  Q.  Snurr  and  Edward  J.  Maginn,  “Molecular  Modeling 
and  Experimental  Studies  of  the  Thermodynamic  and  Transport  Properties  of  Pyridinium- 
Based  Ionic  Liquids”,  /  Physical  Chemistry  B.,  2006,  110,  2821-2832. 

3.2.1  Motivation 

Ionic  liquids  having  a  pyridinium-based  cation  are  another  import  class  of  compounds 
with  potential  use  to  the  Air  Force.  From  other  work  at  Notre  Dame,  we  have  learned  that 
pyridinium-based  ionic  liquids  have  the  highest  thermal  stability  of  any  ionic  liquid 
examined  to  date.  This  property  could  lead  to  many  interesting  applications  for  the  Air 
Force.  As  is  the  case  with  the  triazo I i uni-based  ionic  liquids,  there  are  no  published 
forcefields  for  this  class  of  compounds,  nor  have  any  simulations  been  carried  out  on 
these  compounds.  Therefore,  we  followed  a  procedure  similar  to  that  outlined  in  Section 

3.1  and  developed  a  foreefield  for  pyridinium-based  ionic  liquids.  We  compared 
thermodynamic  properties  against  available  experimental  data.  We  were  also  interested  in 
determining  how  well  lime-dependent  properties  could  be  modeled.  We  therefore  enlisted 
the  help  of  researchers  at  Northwestern  University,  who  used  pulsed-field  gradient 
nuclear  magnetic  resonance  (PFG  NMR)  to  measure  self-diffusivities  for  three  ionic 
liquid  systems.  We  then  compared  our  calculated  self-diffusivities  against  these  data. 


3.2.2  Summary  of  Major  Findings 

The  following  ionic  liquids  were  simulated:  l-n-hexyl-3-methylpyridinium 
bis(trinuoromethanesulfonyi)amide  ([hmpy][Tf:N]),  l-«-hexyl-3,5-dimethylpyridinium 
bis(trifluoromethanesulfonyl)amide  ([hdmpy][Tf;iN]),  and  1 -n-octyl-3-methyl-pyridinium 
bis(trifluoromcthanesulfony[)amide  [ompy][Tf2N], 


Thermodynamic  Properties 

The  original  foreefield  we  developed  will  be  designated  FF1.  Figure  5  shows  the 
computed  density  as  a  function  of  temperature  for  [ompy][Tf:N]  compared  with 
experimental  data  from  our  lab.  As  can  be  seen,  FF1  yields  densities  that  are  abot  5%  too 
high.  To  improve  upon  this,  the  collision  diameters  of  the  ring  atoms  (N,  C  and  H)  were 
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increased  by  15%.  This  set  of  parameters,  denoted  FF2,  does  a  much  better  job  matching 
experimental  densities.  Similar  results  were  obtained  for  the  other  ionic  liquids. 


Figure  5;  Computed  density  as  a  function  of  temperature  for  [ompy][Tf}N],  FF1  is  the 
original  forcefield,  while  FF2  is  the  re-parameterized  forcefteld.  The  line  is  experimental 
data . 


Table  5  lists  computed  and  (when  available)  experimental  values  for  the  volumetric 
expansivity,  isothermal  compressibility,  and  heat  capacities.  The  computed  volumetric 
expansivities  for  both  FF1  and  FF2  are  slightly  lower  than  the  experimental  values, 
ranging  from  4-5x10  4  K  1  while  experimental  values  range  from  6-7x10 4  K'1.  Isothermal 
compressibilities  computed  with  FF1  are  just  under  2  x  10'5  bar1  while  experimental 
values  for  [hmpy][Tf2N]  and  [hdmpy]Tf2N]  are  found  to  be  on  the  order  of  3  x  10'5  bar1. 
This  is  similar  to  the  isothermal  compressibility  of  various  imidazolium-based  ionic 
liquids.  The  general  conclusion  that  may  be  drawn  from  these  results  is  that  the  three 
systems  show  the  same  tendency  to  expand  /  contract  under  the  influence  of  temperature 
and  pressure,  and  thus  the  intermolecular  forces  among  the  various  ions  are  similar.  This 
is  not  surprising,  given  the  structural  similarities  among  the  three  liquids  as  well  as  the 
prominent  role  electrostatics  plays  in  the  volumetric  properties  of  these  liquids. 

Table  5:  Computed  expansivities,  heat  capacities  and  compressibilities  for  pyridinium- 
based  ionic  liquids. 


IL 

(WK-’jxlO4 

atK'VlO4 

a'^K’^xlO4 

[hmpy][Tf2Nj 

6.35i 

5.2338 

3.99„ 

[ompylfTfjNl 

6.942 

4.9539 

4.51 5S 

[hdmpy][Tf2N] 

6.742 

4.02ap 

3.8546 

CP  (J/K  mol) 

Cv  (J/K  mol) 

KT(bar'1)x10l 

[hmpy][Tf2N] 

1291,6 

113533 

1  8421 

[ompy][Tf2N] 

145423 

131036 

1.9722 

[hdmpy][Tf2N] 

1345,3 

125241 

1  -9323 

*  FFI 
**  FF2 
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Heat  capacities  at  constant  volume  and  pressure  range  from  about  1.2  kJ/mol  K.  for 
[hmpyJfTfjN]  to  1.4  kJ/mol  K  for  [ompy][Tf2Nj.  Based  on  results  from  our  lab,  this 
appears  to  be  high  by  about  a  factor  of  two.  Experimental  values  are  closer  to  the  range  of 
0.6-0. 7  kJ/mol  K.  It  is  encouraging,  however,  that  the  forcefield  captures  the 
experimental  trend  that  [ompy][Tf2N]  has  the  highest  heat  capacity,  followed  by 
[hdmpy][Tf2N],  with  [hmpy][Tf2N]  having  the  lowest  heat  capacity.  The  longer  the 
length  of  the  alkyl  chain  or  more  substitution  along  the  ring,  the  more  energy  storage 
mechanisms  the  cation  has,  and  thus  the  higher  the  heat  capacity. 

Dynamic  Properties 

The  microscopic  dynamics  of  ionic  liquids  plays  a  critical  role  in  determining  the 
rheological  properties  of  these  compounds.  It  also  governs  mass  transfer  behavior,  which 
is  crucial  in  determining  chemical  reactivity  in  ionic  liquid  solvents  and  the  performance 
of  ionic  liquids  in  separations  and  electrochemical  applications.  A  particularly  useful 
measure  of  liquid  dynamics  that  is  amenable  to  both  experimental  and  computational 
determination  is  the  self-diffusivity,  Ds,  defined  as 


n  ]v  d 
D  =  -lim  — 

1  6  dt 


(2) 


where  the  term  in  pointed  brackets  is  the  mean-squared  displacement  (MSD)  of  a  tagged 
particle.  Previous  PFG  NMR  studies  of  imidazolium-based  ionic  liquids  found  that  the 
cation  has  a  slightly  higher  self-diffusivity  than  the  anion,  despite  its  larger  size8.  These 
studies  also  found  that  the  self — diffusivity  correlates  well  with  the  inverse  of  the 
viscosity,  but  application  of  the  Stokes-Einstein  model  yields  ion  sizes  that  are  not 
qualitatively  correct. 

Our  collaborator  at  Northwestern  University  (Prof.  Randall  Snun)  performed  the  PFG 
NMR  measurements  and  we  computed  the  self-diffusivity  using  MD.  Table  6  lists  the 
results.  Two  main  observations  can  be  made  from  the  experimental  data. 

First,  the  experimental  self-diffusivitics  of  the  cations  and  anions  are  similar  for  each 
ionic  liquid,  and  are  on  the  order  of  1x10'"  m2/s  at  room  temperature.  This  is  roughly  two 
orders  of  magnitude  smaller  than  that  of  water.  [hmpy][Tf2N]  has  the  highest  self- 
diffusivity,  while  [ompy][TfzN]  has  the  lowest,  with  [hdmpy][Tf2N]  falling  between 
these  two.  This  trend  is  consistent  with  the  molar  volumes  and  viscosities;  the  higher  the 
molar  volume  and  lower  the  viscosity  of  the  liquid,  the  higher  the  self-diffusivity. 

Second,  for  [hmpy][Tf2N]  the  cation  appears  to  have  a  slightly  higher  self-diffusivity 
than  the  anion  at  all  temperatures,  although  given  the  uncertainty  in  the  data  this 
difference  is  not  statistically  significant.  For  the  other  two  ionic  liquids,  the  anion  clearly 
has  the  larger  self-diffusivity.  The  bulky  [ompy]  cation  has  the  lowest  self-diffusivity  at 
all  temperatures,  while  [hdmpy]  has  a  self-diffusivity  mid-way  between  [ompy]  and 
[hmpy].  As  noted  above,  previous  PFG  NMR  studies  of  imidazolium-based  ionic  liquids 
always  found  that  the  cations  had  larger  self-diffusivities  than  their  smaller  anion  pairs. 
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Table  6:  Experimental  and  computed  self-diffusivities  as  a  function  of  temperature  for  the 
cations  and  anions. 


PFG  NMR  Self  Difusivitiesx10"(m2/s) 

Of*  1  Ll 


ec 

,bF 

<  exp 

He*P  (Avrg) 

Anionsim 

Cationsim 

[hmpy][Tf;N] 

25 

1.207 

1.21 8005 

0.083 

0.108 

35 

2.122 

2,145oi3 

0.174 

0.197 

45 

3.415 

3.4480i7 

0.281 

0.329 

55 

5.544 

5.590042 

0.214 

0.278 

65 

9.658 

9.61  3oab 

0.445 

0.492 

75 

17.690 

17.867,97 

0.626 

0.672 

100 

1.548 

1.448 

125 

3.206 

3.651 

150 

6.219 

6.737 

[ompy][TfzN] 

25 

0.916 

O.86I004 

0.150 

0,144 

35 

1.659 

1  .555o23 

0.172 

0.186 

45 

2.713 

2.530025 

0.216 

0.255 

55 

4.302 

4.01  1q23 

0.343 

0.325 

65 

7.158 

6.731 037 

0.423 

0.558 

75 

12.400 

1 1 7Q7i2i 

0,616 

0.631 

100 

1.221 

1.269 

125 

2,620 

3.248 

150 

5.657 

5.801 

[hdmpy][Tfz  NJ 

25 

0.980 

0.925007 

0.127 

0.119 

35 

1.804 

1  ♦674oos 

0.196 

0.202 

45 

3.017 

2.844™ 

0.185 

0.226 

55 

5,163 

4.867o76 

0.230 

0.268 

65 

8.847 

8 .450,  ,9 

0.380 

0.378 

75 

17.130 

1  7.090268 

0.506 

0.589 

100 

1.384 

1.380 

125 

2.596 

2.697 

150 

6.871 

6.597 

The  present  results  suggest  that  this  is  not  a  universal  phenomenon,  but  that  the  nature  of 
the  ions  themselves  is  important  in  determining  the  relationship  between  cation  and  anion 
self-diffusivities.  Increasing  alkyl  substitution  and  chain  length  apparently  serve  to 
hinder  cation  diffusion  in  these  systems  such  that  the  anions  now  diffuse  faster  than  the 
cations.  It  is  interesting  to  note,  however,  that  anion  diffusivities  do  track  cation 
diffusivities.  That  is,  the  [TfzNJ  anion  has  the  highest  sclf-diffusivity  when  paired  with 
the  fastest  diffusing  [hmpy]  cation,  and  the  lowest  self-diffusivity  when  paired  with  the 
slowest  diffusing  [ompy]  cation.  This  suggests  that  there  is  ion  pairing  present  in  the 
liquid  such  that  ions  do  not  exhibit  free  diffusion,  but  rather  travel  in  pairs  or  clusters. 
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To  gain  additional  insight  into  the  operative  diffusion  mechanisms,  the  self-diffusivity 
was  computed  by  applying  eqn  2  to  the  results  of  microcanonical  ensemble  molecular 
dynamics  simulations.  It  is  important  to  note  that  eqn  2  is  valid  only  when  true  diffusive 
motion  is  observed.  For  the  PFG  NMR  experiments,  the  diffusion  times  over  which  data 
were  taken  were  very  long  (63  ms),  so  that  diffusional  motion  was  assured.  Such  long 
times  are  not  attainable  with  molecular  dynamics,  however,  so  great  care  must  be  taken  to 
ensure  that  true  diffusive  processes  are  simulated.  Fig.  6  shows  the  MSD  as  a  function  of 
time  for  the  [ompy]  cation  in  [ompy][Tf2N]  at  298  K  and  423  K.  The  solid  line  is  the  total 
MSD  while  the  dashed  lines  are  for  displacements  along  the  three  Cartesian  directions. 
These  latter  values  give  an  indication  of  the  statistical  variation  in  the  MSD  and  also 
confirm  that  there  is  no  preferential  diffusion  along  a  given  Cartesian  axis.  Both  plots 
appear  to  show  linear  behavior  at  long  times,  and  based  on  this  one  would  be  tempted  to 
directly  apply  eqn  2  to  compute  Ds.  Upon  closer  inspection  of  the  low  temperature 
results,  however,  one  finds  that  this  system  is  not  yet  in  the  diffusive  regime. 


Figure  6:  MSD  for  [ompy]  cation  as  a  function  of  time  at  298  K  (left)  and  423  K  (right). 


The  dynamic  nature  of  a  liquid  can  be  characterized  by  examining  the  way  in  which  the 
mean-squared  displacement  scales  with  time.  This  can  be  quantified  as 


«  p 


(3) 


where  P  characterizes  the  type  of  motion  present  in  the  system.  At  very  short  times  when 
ballistic  motion  dominates,  p=2.  For  long-time  diffusive  motion,  eqn  2  implies  (1=1.  At 
intermediate  time  scales,  however,  sub-diffusive  dynamics  characteristic  of  glass-like 
behavior  can  result  in  P  <  I,  a  situation  observed  previously  with  ionic  liquids4.  By 
plotting  the  following  expression  as  a  function  of  time 


m= 


rflog(Ar2} 

d\oz(t) 


(4) 
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one  can  easily  determine  which  dynamic  regime  a  system  is  in.  Fig.  7  shows  that  at  298 
K,  the  system  has  not  yet  reached  diffusive  behavior  even  after  5  ns;  fl  has  only  reached  a 
value  of  0.6.  At  the  higher  temperature,  the  system  does  approach  diffusive  behavior, 
with  (J  exceeding  0.9  after  4.5  ns.  Note  that  the  resulting  self-diffusivities  depend  strongly 
upon  the  time  interval  over  which  the  MSD  is  taken. 


Figure  7:  Slope  parameter  as  a  function  of  time  for  [ompyj  cation  in  [ompy][Tf2N]  at 
298  K  (left)  and  423  K  (right).  Diffusive  motion  is  acheived  when  beta  =  J. 


These  results,  coupled  with  other  analyses  we  have  performed,  indicate  that  these  systems 
exhibit  glass-like  dynamics  at  room  temperature,  even  though  they  are  technically  above 
their  glass  transition  temperature.  Because  of  this,  computing  self-diffusivities  at  low 
temperatures  requires  extraordinarily  long  simulations.  We  have  extended  these 
calculations  for  20  ns  and  still  do  not  see  linear  MSD  behavior  at  298  K.  This  suggests 
that  we  must  be  careful  when  reporting  self — diffusi vites  at  low  temperature. 

Clearly,  the  forcefield  we  developed  matches  the  thermodynamic  properties  of  these 
systems  reasonably  well,  but  seriously  underestimates  the  mobility.  Similar  behavior  has 
been  observed  in  simulations  of  alkali  halide  salts.  Part  of  this  may  be  due  to  a  neglect  of 
polarizabiltiy  in  the  forcefield,  as  suggested  by  Voth  and  co-workers9.  Further  study  of 
this  issue  is  warranted. 

3.2  J  Relationship  to  the  Goals  of  the  Project. 

As  was  the  case  with  the  triazolium  study,  the  development  of  forcefields  that  can  be 
used  to  accurately  predict  thermodynamic  and  transport  properties  of  ionic  liquids  is 
essential.  We  have  developed  a  new  forcefield  for  pyridinium-based  ionic  liquids,  which 
show  excellent  potential  for  Air  Force  applications  due  to  their  high  thermal  stability.  We 
also  have  much  more  experimental  data  for  these  systems,  which  enables  careful 
benchmarks  to  be  made.  We  have  shown  that  we  can  make  accurate  predictions  of 
volumetric  properties,  but  timc-dcpcndcnt  properties  (such  as  diffusivity)  arc  more 
difficult  and  will  require  additional  effort. 
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3.3  Development  of  New  Simulations  Methods  for  Computing 
Solubility  in  Ionic  Liquids 

Jindal  K.  Shah  and  Edward  J.  Maginn,  “A  Monte  Carlo  Simulation  Study  of  the  Ionic 
Liquid  l-n-Butyl-3-Methylimidazolium  Hexafluorophosphate:  Liquid  Structure, 
Volumetric  Properties,  and  Infinite  Dilution  Solution  Thermodynamics  of  CO2",  Fluid 
Phase  Equilibria,  2004,222-223,  195-203, 

Cesar  Cadcna,  Jennifer  L.  Anthony,  Jindal  K.  Shah,  Timothy  1.  Morrow,  Joan  F. 
Brennecke  and  Edward  J.  Maginn,  “Why  is  CO2  So  Soluble  in  Imidazolium-based  Ionic 
Liquids?”,  Journal  of  the  American  Chemical  Society,  2004,  126,  5300-5308. 

Jindal  K.  Shah  and  Edward  J,  Maginn,  “Monte  Carlo  Simulations  of  Gas  Solubility  in  the 
Ionic  Liquid  I-«-butyl-3-methylimidazolium  hexafluorophosphate”,  J.  Physical 
Chemistry  B.,  2005,  109,  10395-10405. 

T.  I.  Morrow  and  E.  J.  Maginn,  "Isomolar-semi grand  ensemble  molecular  dynamics: 
Application  to  vapor-liquid  equilibrium  of  the  mixture  methane/ethane".  Journal  of 
Chemical  Physics,  2006,  125,  204712. 


3.3.1  Motivation 

Any  practical  use  of  an  ionic  liquid  will  involve  exposure  to  the  atmosphere.  It  is  known 
from  experiment  that  ionic  liquids  tend  to  dissolve  water,  CO2,  and  other  species  from  the 
atmosphere.  In  fact,  CO2  solubilities  are  so  high  that  ionic  liquids  are  being  investigated 
for  CO2  capture  applications.  When  gases  dissolve  in  ionic  liquids,  the  properties  of  the 
ionic  liquids  can  change  dramatically.  Hence,  it  is  important  to  not  only  understand  the 
properties  of  pure  ionic  liquids,  but  it  is  equally  important  to  examine  properties  of 
mixtures.  We  would  also  like  to  understand  the  factors  that  govern  whether  an  ionic 
liquid  has  a  high  or  low  affinity  for  a  particular  species.  Because  of  this,  we  were 
motivated  to  carry  out  calculations  of  the  solubility  of  various  gases  in  ionic  liquids. 
Standard  simulation  techniques  for  doing  this  turned  out  to  be  inadequate,  so  we  have 
worked  on  the  development  of  new  simulation  methods  for  carrying  out  these 
calculations. 

3.3.2  Summary  of  Major  Findings 

Simulating  Henry’s  Law  Constant  of  Gases  in  Ionic  Liquids 

The  Henry’s  Law  constant,  H,  is  defined  as 

//slim-  (5) 

p-,0  x 

where /  is  the  fugacity  of  the  gas  and  x  is  the  mol  fraction  of  the  gas  in  the  liquid.  The 
Henry’s  Law  constant  thus  describes  the  low  solubility  limit,  and  is  directly  related  to  the 
free  energy  of  solvation.  Several  techniques  may  be  used  to  compute  H.  One  way  is  to 
compute  the  isotherm  of  the  species  of  interest,  and  then  take  the  limiting  slope.  Another 
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way  is  to  compute  the  free  energy  of  solvation  and  relate  this  to  H.  We  have  used  both 
methods,  as  summarized  below. 

The  Henry's  Law  constant  is  related  to  the  excess  chemical  potential  fi* *  via  the 
following  relation 


H  =  pRTexp[Pne’] 


(6) 


where  pis  the  liquid  density  and  pis  inverse  temperature.  The  excess  chemical  potential 
is  traditionally  computed  using  Widom’s  “test  particle  insertion11  method,  in  which  “test” 
gas  molecules  are  randomly  inserted  into  the  fluid  and  the  ensemble  average  energy 
tallied.  In  the  isothermaLisobaric  ensemble,  the  excess  chemical  potential  is  obtained  by 
the  following  expression 


p"  =  —kT 


(Vexp  [-ffl/.JL 


(V) 


SFT 


(7) 


where  Ug  is  the  energy  of  the  test  or  “ghost”  particle.  We  used  this  approach  and,  as 
shown  below,  (he  method  is  not  always  reliable  due  to  the  strong  overlap  that  is  present 
between  the  gas  solute  and  the  dense  liquid.  For  that  reason,  we  also  developed  a  new 
“expanded  ensemble”  (EE)  approach  for  computing  H. 

The  EE  method  essentially  overcomes  the  problems  associated  with  the  test  particle 
method  by  gradually  “growing”  the  gas  molecule  into  the  system  instead  of  attempting  to 
insert  the  molecule  all  at  once.  The  intermolecular  potential  energy  between  the  gas 
molecule  and  the  ionic  liquid  is  gradually  turned  on  by  stochastically  modifying  a 
coupling  parameter.  By  recording  the  fraction  of  the  time  a  solute  spends  in  the  fully 
coupled  state  (A  -  l)and  a  fully  decoupled  state  (A  =  0),  one  can  directly  determine  the 
free  energy  of  solvation  and  hence  the  Henry’s  Law  constant. 
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Figure  8:  Convergence  of  the  Henry's  Law  constant  for  oxygen  in  [bmim][PF6]  using  the 
test  particle  insertion  method. 

Figure  8  shows  the  results  of  a  test  particle  insertion  Monte  Carlo  calculation  at  298  K 
and  1  bar.  It  demonstrates  how  the  test  particle  insertion  method  converges  for  O2  in  1  -«■ 
butyl-3-methylimidazolium  hexafluorophosphate  ([bmim][PFe]).  Each  individual 
trajectory  appears  to  converge,  but  the  variation  between  any  given  simulation  is  fairly 
high.  In  Figure  9,  the  same  result  is  shown  for  the  EE  method.  Although  the  convergence 
trends  look  similar  to  that  with  the  particle  insertion  method,  we  find  that  the  EE  method 
yields  more  reliable  results  and  does  not  depend  as  strongly  on  system  size  as  the  test 
particle  method.  For  larger  gas  molecules  such  as  ethane,  the  advantage  of  the  EE  method 
is  even  grerater.  Another  advantage  of  the  EE  method  is  that,  since  the  solute  molecules 
actually  interact  with  the  fluid,  structural  features  of  the  solute  dissolved  in  the  solvent 
can  be  obtained. 
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Figure  9:  Convergence  of  EE  simulations  for  O:  in  [bmim][PF6]. 


Figure  10  shows  site-site  radial  distribution  functions  for  the  0  and  H  atoms  in  water 
with  the  fluorine  and  phosphorous  atoms  of  the  anion  and  the  acidic  C2  carbon  of  the 
imidazolium  ring.  The  peak  in  the  F-H  trace  at  0,18  nm  is  a  clear  signature  of  hydrogen 
bonding  between  the  anion  and  water.  The  association  with  the  carbon  atom  of  the  cation 
is  much  weaker.  This  suggest  that  water  can  H-bond  with  the  anion,  and  that  the  nature  of 
the  anion  should  dominate  water  solubility.  In  fact,  this  is  exactly  what  is  seen 
experimentally.  By  changing  the  nature  of  the  anion,  one  can  control  the  water  solubility 
to  a  much  larger  extent  than  is  possible  by  changing  the  cation. 


Figure  10:  Radial  distribution  functions  for  water  associating  with  various  sites  on  the 
cation  and  anion  of  [ bmimJfPFi ].  Clear  evidence  of  hydrogen  bonding  is  observed. 
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Development  of  a  Semi-Grand  Molecular  Dynamics  Methods  for  Computing 
Solubilities  in  Liquids 

An  alternative  to  EE  and  test  particle  insertion  methods  is  to  use  a  so-called  semi-grand 
ensemble  technique.  In  this  method,  the  identity  of  a  molecule  is  changed  from  one 
species  to  another,  typically  using  using  a  stochastic  procedure.  This  only  really  works 
when  the  size  and  nature  of  the  two  species  are  similar.  When  the  size  and  /  or  energetics 
of  the  species  in  the  mixture  are  very  different  (as  is  the  case  of  small  solutes  in  ionic 
liquids)  then  this  method  fails  because  random  changes  in  the  “identity”  of  a  molecule 
will  usually  lead  to  molecules  overlapping  each  other  in  high-energy  configurations.  To 
overcome  this  difficulty,  we  recently  developed  a  dynamic  method  for  exchanging  the 
identity  of  a  molecule".  With  the  goal  of  applying  this  technique  to  ionic  liquid 
simulations,  we  first  implemented  it  in  the  simplest  possible  case:  a  mixture  of  methane 
and  ethane.  We  call  this  new  technique  isomolar  semi-grand  ensemble  molecular 
dynamics  (iSGMD). 

In  iSGMD  the  transformation  of  molecules  occurs  gradually  and  dynamically  by  adding 
an  extended  system  variable  to  the  normal  isothermal-isobaric  equations  of  motion.  The 
role  of  the  new  extended  system  variable  is  to  dynamically  transform  a  molecule  between 
two  identities.  The  dynamical  transformation  scheme  allows  the  system  to  automatically 
adjust  to  the  difference  between  the  setpoint  and  instantaneous  chemical  potential 
differences  between  any  species  i  and  reference  species  \,  (ji, ,  -  jii),  under  the  influence 
of  the  forces  present  in  the  system.  The  equations  of  motion  for  carrying  this  out  arc  quite 
complex,  and  can  be  found  in  our  publication".  Instead  of  specifying  a  chemical 
potential  difference,  it  is  more  convenient  to  specify  a  fugacity  fraction,  defined  as 

(,  -  (8) 

X/, 

i 

In  a  binary  mixture,  the  fugacity  fraction  £wi]|  range  from  0-1.  Once  this  value  is 
specified,  the  iSGMD  dynamical  equations  will  begin  transforming  the  identity  of  the 
molecules  until  a  steady-state  composition  is  reached.  By  knowing  how  composition 
varies  with  fugacity  fraction,  one  can  use  Gibbs- Duhem  integration  methods  to  map  out 
an  entire  phase  diagram  or  isotherm.  Figure  1 1  shows  the  results  for  the  simple  test 
system  of  methane-ethane.  As  can  be  seen,  the  method  yields  results  that  are  in  excellent 
agreement  with  experiment12  and  with  a  recent  Monte  Carlo13  simulation  of  the  same 
system.  Importantly,  the  swapping  efficiency  of  our  new  method  is  a  factor  of  ten  greater 
than  the  standard  Monte  Carlo  approach.  While  the  inefficiency  of  the  conventional 
Monte  Carlo  approach  does  not  prevent  its  use  for  this  simple  system,  it  will  not  be 
possible  to  use  that  method  on  ionic  liquids.  The  improved  efficiency  of  our  new  iSGMD 
method,  however,  means  that  we  can  now  use  it  to  begin  tackling  phase  equilibria 
problems  with  ionic  liquids.  We  hope  to  extend  this  method  to  the  investigation  of  ionic 
liquids  in  the  near  future. 
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Figure  11:  Computed  pressure  versus  fugacity  fraction  diagram  for  methane  -  ethane. 
Experimental  data  shown  as  the  line,  solid  symbols  are  our  results,  and  open  symbols  are 
previous  MC  results. 

Development  of  a  Semi-Grand  Ensemble  Hybrid  Monte  Carlo  Method  for 
Computing  Gas  Absorption  Isotherms 

Finally,  we  describe  a  new  method  we  have  developed  which  we  believe  will  be  the 
easiest  to  apply  for  computing  absorption  isotherms  of  small  molecules  in  ionic  liquids. 
This  method  utilizes  a  semi-grand  ensemble  in  which  the  total  number  of  ionic  liquid 
molecules  is  fixed,  as  is  the  total  pressure,  temperature  and  the  fugacity  (partial  pressure) 
of  the  gas  species.  The  concentration  of  the  gas  solute  in  the  liquid  phase  adapts  to  the 
imposed  fugacity,  thus  yielding  a  point  on  an  absorption  isotherm. 

Equilibration  of  the  liquid  phase  is  accomplished  via  a  hybrid  Monte  Carlo  method  in 
which  short  microcanonical  ensemble  MD  trajectories  are  run,  and  the  Hamiltonian, 
defined  as 


N 


ffM=Z5r+£/M 


/=!  2 m. 


(9) 


is  computed  before  and  after  the  trajectory.  An  equilibration  “move”  is  then  accepted 
according  to  the  following  rule 


P  =min(l,e-M//). 

QCCfrans  \  ’  / 


(10) 


Volume  moves  are  required  to  maintain  a  constant  pressure.  These  are  done  in  a 
Metropolis-like  fashion,  with  acceptance  probability  given  by 
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Gas  molecules  are  inserted  and  deleted  into  the  system  in  a  biased  manner,  with  the 
acceptance  probability  for  insertion  given  by 


P  ,  =  min 

acc,  insert 


i. 


a 

_ n 
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&  f(5V 


Z'g  N 


-exp[-p(un~um)] 


(12) 


and  deletions 
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(13) 


where  a  is  a  cavity-biasing  factor,  the  form  of  which  is  given  elsewhere14.  The  cavity¬ 
biasing  algorithm  is  needed  to  improve  the  acceptance  probabilities  of  insertion.  The 
method  was  tested  by  reproducing  previous  simulation  results  for  argon  and  water  as  well 
as  neat  [bmimJJPFe].  The  method  was  then  applied  to  compute  the  absorption  isotherm  of 
CO2  in  [bmimjfPFe],  The  initial  result  of  this  calculation  looks  very  promising;  Figure  12 
shows  the  computed  points  on  the  isotherm,  along  with  experimental  data  of  Akils.  So 
far,  the  results  are  very  encouraging  and  suggest  that  we  can  reliably  compute  gas 
solubilities  in  ionic  liquids. 


160 

140 

*  Aki  et  al. 

120  *  Simulated 


ra 

n 

2? 

□  BO 


8 

Cl  60 


0 

0  0,1  0.2  0.3  04  0.5  06  0.7 

C03t  mole  fraction 

Figure  12:  Experimental  data  (Aki)  and  computed  isotherm  for  CO2  in  [bmim][PF6]  at 
298  K 
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3.3.3  Relationship  to  the  Coals  of  the  Project 

We  have  made  extensive  progress  in  developing  methods  that  can  compute  the 
equilibrium  composition  of  solutes  dissolved  in  ionic  liquids.  This  ability  is  crucial,  since 
in  actual  applications,  ionic  liquids  will  be  exposed  to  water  and  permanent  atmospheric 
gases.  How  much  of  these  species  will  dissolve  in  an  ionic  liquid  is  difficult  to  predict. 
We  do  know  that  the  properties  of  the  liquid  can  be  greatly  altered  by  these  dissolved 
species.  Therefore,  the  methods  developed  here  are  of  critical  importance  to  our  efforts  to 
learn  how  to  model  these  systems.  Given  these  tools,  we  can  compute  the  equilibrium 
composition  of  species  dissolved  in  an  ionic  liquid  under  application  conditions.  With 
this  information  and  the  forcefields  we  have  developed  to  date,  we  can  then  perform 
mixture  simulations  to  assess  how  the  properties  of  the  liquid  will  be  influenced  by  the 
solutes.  An  example  of  how  this  can  be  done  is  given  in  the  next  section.  These  powerful 
new  methods  are  not  limited  to  ionic  liquids  calculations,  however,  but  can  be  used  on  a 
wide  range  of  complex  fluid  systems. 


3.4  Simulating  mixtures  containing  imidazolium-  and  pyridinium- 
based  ionic  liquids  and  1 -butanol 

J.  M.  Crosthwaite,  S.  N.  V.  K.  Aki,  E.  J.  Maginn  and  J.  F.  Brcnnecke,  “Liquid  Phase 
Behavior  of  Imidazolium-Based  Ionic  Liquids  with  Alcohols:  Effect  of  Hydrogen 
Bonding  and  Non-Polar  Interactions”,  Fluid  Phase  Equilibria,  2005,  228,  303-309. 

T.  I.  Morrow  and  E.  J.  Maginn,  “Molecular  Simulation  of  Mixtures  Containing 
Imidazolium  and  Pyridinium-based  Ionic  Liquids  and  1 -Butanol”,  ACS  Symposium 
Series,  in  press. 

3.4.1  Motivation 

It  is  known  from  the  experiments  of  our  collaborator  Prof.  Joan  Brennecke  that  many 
ionic  liquid  /  alcohol  mixtures  form  two  immiscible  liquid  phases.  In  particular, 
imidazolium-  and  pyridinium-based  ionic  liquids  have  been  found  to  phase  separate  into 
two  liquid  phases  when  mixed  with  1 -butanol.  Given  that  we  have  forcefields  for  all 
these  compounds,  we  decided  to  carry  out  simulations  to  see  if  we  could  understand  this 
behavior  better. 

3.4.2  Summary  of  Major  Findings 

Solutions  of  I -butanol  and  l-n-butyl-3-methyIimidazolium  tetrafluoroborate 
([bmim][BF4]),  l-n-butyl-3-methylimidazolium  bis(trifluoromethanesulfonyl)  amide 
([bmim][Tf2N]),  l-n-butyl-3-methylpyridinium  tetrafluoroborate  ([bmpyjfBF^]),  and  1-n- 
butyl-3-methylpyridinium  bis(trifluoromethanesulfonyl)  amide  ([bmpy][Tf2N])  were 
examined  using  isothermal-isobaric  molecular  dynamics  simulations.  Quantities 
computed  include  molar  volumes,  self-difftisivities,  radial  distribution  functions,  local 
composition  functions,  and  interaction  energies.  We  were  particularly  interested  in  using 
the  simulations  to  help  explain  recent  experimental  results  in  which  [bmpyjiBF^]  was 

found  to  have  a  lower  upper  critical  solution  temperature  (UCST)  with  1 -butanol  than 
[bmimjfBF^.  In  addition,  it  was  found  that  [bmim][Tf2N]  has  a  lower  UCST  with  I- 
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butanol  than  [bmpy][Tf2N].  Apparently,  the  individual  nature  of  the  cation  and  anion  is 
not  enough  to  explain  such  trends,  but  one  must  consider  all  the  interactions  present. 


Simulations  at  two  compositions  were  performed  for  each  ionic  liquid  /  butanol  mixture, 
with  one  simulation  consisting  of  a  butanol-rich  system  (-4  mol%  ionic  liquid),  and  the 
other  simulation  consisting  of  an  IL-rich  system  (~25-40  mol%  ionic  liquid).  These 
compositions  are  all  near  the  experimental  coexistence  compositions.  All  eight 
simulations  were  performed  at  P=1  bar  and  a  reduced  temperature. 


=  0.9836,  where  Tucst  is  the  experimentally  determined  UCST. 


We  computed  thermodynamic  and  transport  quantities  for  the  mixture,  as  well  as  the 
liquid  structure  and  energetic  interactions.  The  net  result  of  this  was  that  we  could  in 
fact  explain  the  trends  in  the  UCST.  Perhaps  the  most  important  factor  is  the  energy  of 
interaction  among  the  various  species  in  the  mixture.  Table  7  shows  the  breakdown  of  the 
potential  energies  of  the  ionic  liquid-rich  phases  into  contributions  from  butanol-butanol, 
butanol-ion,  and  ion-ion  interactions  for  [bmim][BF4]  and  [bmpy][BF4].  The  butanol- 
butanol  interaction  is  1.0  kcal  /  mol  more  favorable  in  the  [bmpy]  mixture  than  the 
[bmim]  mixture.  On  the  other  hand,  the  ion-ion  interactions  are  12.7  kcal  /  mol  more 
favorable  in  the  [bmim]  mixture,  which  demonstrates  that  the  [bmim][BF4]  interaction  is 
stronger  than  that  of  the  [bmpy][BF4]  interaction.  The  butanol-ion  interactions  are  only 
slightly  (0.4  kcal  /  mol)  more  favorable  in  the  [bmpy]  mixture  than  in  the  [bmim] 
mixture.  Altogether  butanol  experiences  about  1.4  kcal  /  mol  lower  energy  in  the 
[bmpy][BF4]  mixture  than  in  the  [bmim][BF4]  mixture,  which  explains  why  it  has  a  lower 
UCST  (i.e.  it  is  more  favorably  solvated  by  the  alcohol).  Similar  trends  in  the  energies 
are  observed  for  the  alcohol-rich  phase. 


Table  7:  Comparison  of  the  contribution  of  the  butanol-butanol,  butanol-ion,  and  ion-ion 
interactions  to  the  system  total  potential  energy  for  the  IL-rich  phases  of  the 
[bmim|[BF4]  and  [bmpy][BF4]  mixtures  with  1-butanol. 


[bmim] 

(bmpy] 

Interaction 

Energy 

(kcal  /  mol  BuGH) 

Energy 

(kcal  /  mol  BuOH) 

BuOH- 

BuQH 

0.486 

-0.514 

BuOH  -  IL 

-14.9 

-153 

1L-IL 

-50  8 

-38.1 
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Table  8:  Comparison  of  the  contribution  of  the  butanol-butanol,  butanol-ion,  and  ion-ion 
interactions  to  the  system  total  potential  energy  for  the  IL-rich  phases  of  the 
[bmim][TfN]  and  [bmpy][Tf2N]  mixtures  with  I -butanol. 


[bmim] 

[bmpy] 

Interaction 

Energy 

(kcal  /  mol  BuOH) 

Energy 

(kcal  /  mol  BuOH) 

BuOH- 

BuOH 

-1.21 

-1.56 

BuOH  -  IL 

-15J2 

“13.78 

IL-  IL 

-5.8 

-6.2 

Table  8  shows  the  breakdown  of  the  potential  energies  of  the  IL-rich  phases  into 
contributions  from  butanol-butanol,  butanol-ion,  and  ion-ion  interactions  on  a  per  mol 
butanol  basis  for  the  simulations  with  the  [Tf2N]  anion.  The  butanol-ion  interactions  are 
1.34  kcal  /  mol-butanol  more  favorable  in  the  [bmim]  mixture  than  the  [bmpy]  mixture, 
which  suggests  that  the  cation/butanol  interaction  is  stronger  in  the  [bmim]  mixture.  On 
the  other  hand,  the  butanol-butanol  interaction  is  0.35  kcal  /  mol  more  favorable  in  the 
[bmpy]  mixture  than  the  [bmim]  mixture  and  the  ion-ion  interactions  are  0.5  kcal  /  mol 
more  favorable  in  the  [bmpy]  mixture.  Overall,  butanol  has  -1.0  kcal  /  mol  more 
favorable  interactions  with  [bmim][Tf2N]  than  [bmpy][TfzN]  due  to  the  favorable 
cation/butanol  energetics.  This  is  in  contrast  to  the  case  when  [BF4]  was  the  anion,  and 
explains  why  for  this  system  the  UCST  is  lower  for  [bmim][Tf2N]  than  [bmpy][Tf2N]. 
These  differences  arc  driven  by  the  fact  that  the  more  localized  charge  on  [BF4]  makes 
the  cation/anion  interactions  very  favorable,  thus  inhibiting  association  between  the 
alcohol  and  the  IL.  This  is  why  the  UCST  is  higher  for  each  1L  with  a  [BF4j  anion  than 
for  the  [TfzN]  ILs.  On  the  other  hand,  the  more  delocalized  charge  on  [Tf2N]  weakens 
the  cation/anion  interactions,  enabling  the  alcohol  to  interact  more  directly  with  the  ions. 
In  this  case,  the  butanol  interacts  most  strongly  with  the  localized  charge  on  [bmim]. 
This  tends  to  lower  its  UCST  relative  to  [bmpy][Tf2N]. 

3.4.3  Relationship  to  the  Goals  of  the  Project 

Given  that  we  can  compute  solubilities  of  species  in  different  ionic  liquids,  we  still  need 
to  be  able  to  assess  how  the  presence  of  dissolved  species  affects  the  properties  of  the 
mixture.  In  this  work,  we  have  examined  the  liquid-liquid  equilibria  of  ionic  liquid 
mixtures.  We  have  also  computed  volumetric,  thermal  and  transport  properties  of  the 
mixture  (not  shown  here).  We  can  use  this  same  formalism  to  examine  other  mixtures  as 
appropriate. 
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3.5  Development  of  a  general  and  robust  simulation  method  for 
predicting  melting  points 

David  M.  Eike,  Joan  F.  Brennecke,  and  Edward  J.  Maginin,  “Toward  a  Robust  and 
General  Molecular  Simulation  Method  for  Computing  Solid-Liquid  Coexistence", 
Journal  of  Chemical  Physics,  2005,122,014115-1-014115-12. 

David  M.  Eike  and  Edward  J.  Maginn,  “Atomistic  Simulation  of  Solid-Liquid 
Coexistence  for  Molecular  Systems:  Application  to  Triazole  and  Benzene",  Journal  of 
Chemical  Physics,  2006,  124,  164503. 

3.5.1  Motivation 

The  ability  to  predict  melting  points  for  ionic  liquids  from  purely  molecular-level 
information  is  of  great  fundamental  and  practical  importance.  Our  understanding  of  what 
makes  these  liquids  low-melting  (or  unexpectedly  high  melting)  is  incomplete.  Given  the 
huge  number  of  potential  compounds  that  can  be  made,  it  is  obviously  inefficient  and 
impractical  to  synthesize  every  possible  compound  in  the  hope  that  one  has  a  low  enough 
melting  point  for  use.  Theory  and  simulation  can  greatly  assist  in  identifying  likely 
candidates. 

Current  molecular  simulation  methods  are  not  well-suited  for  computing  melting  points 
of  ionic  liquids.  The  standard  method  for  computing  melting  behavior  is  to  separately 
compute  the  absolute  free  energy  of  the  crystalline  and  liquid  phases  by  performing 
thermodynamic  integration  between  a  reference  state  and  the  state  point  of  interest.  These 
free  energies  are  then  matched  and  the  point  at  which  the  solid  and  liquid  free  energies 
are  equal  is  by  definition  the  melting  point.  There  are  a  number  of  significant  technical 
hurdles  that  must  be  overcome  to  apply  this  method  to  ionic  liquids.  We  believe  that 
reference  state  methods  are  not  the  best  way  to  approach  the  problem. 

Another  method  that  has  been  used  successfully  is  to  gradually  heat  a  crystal  and  watch 
for  signatures  of  a  first  order  phase  transition.  For  example,  in  the  isothermal-isobaric 
ensemble,  the  density  will  show  a  sharp  transition  at  the  melting  point.  It  is  observed, 
however,  that  the  melting  point  is  always  overpredictcd  with  this  approach,  due  to  the 
finite  free  energy  barrier  that  must  be  overcome  to  nucleate  a  liquid  phase.  Thus,  this 
method  docs  not  measure  the  true  thermodynamic  melting  point.  To  overcome  this, 
Thompson’s  group  has  shown  that  the  introduction  of  void  defects  into  the  crystal 
reduces  the  melting  point.  At  some  defect  density,  the  simulated  melting  point 
approaches  the  experimental  melting  point.  While  this  approach  is  easy  to  implement  and 
can  give  good  results,  there  are  cases  where  no  clear  melting  point  is  observed. 

To  overcome  these  difficulties,  we  decided  to  develop  a  rigorous  free  energy-based 
simulation  approach  that  can  directly  determine  melting  points  for  arbitrarily  complex 
molecules.  A  brief  summary  of  the  method  and  results  is  given  below. 
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3.5.2  Summary  of  Major  Findings 

The  procedure  for  computing  the  solid— liquid  coexistence  curve  consists  of  five  steps: 
I )  Locate  the  approximate  melting  point  (if  unknown)  using  a  gradual  heating  method;  2) 
Compute  relative  equations  of  state  for  the  solid  and  liquid  phases  using  a  multiple 
histogram  re-weighting  procedure;  3)  Determine  the  difference  in  free  energy  between 
the  solid  and  liquid  phases  at  a  single  statepoint  using  a  newly  developed  pseudo- 
supercritical  inter-phase  path  sampling  method;  4)  Find  the  coexistence  point  by  using 
the  result  from  step  3  to  reference  the  equations  of  state  computed  in  step  2  to  each  other; 
5)  Given  a  single  coexistence  point,  compute  the  coexistence  curve  over  a  range  of  state 
points  using  Gibbs— Duhem  integration. 

This  method  was  tested  by  computing  the  melting  point  for  a  Lennard-Jones  system  as 
well  as  NaCI.  A  brief  summary  of  some  of  these  results  is  given  here.  A  detailed 
description  of  the  method  has  been  submitted  for  publication,  and  can  be  obtained  from 
the  PI  upon  request. 


The  approximate  melting  point  for  sodium  chloride  was  computed  using  the  gradual 
heating  method.  Fig.  13  shows  the  results  of  this  calculation.  It  can  be  seen  that  the 
calculations  overpredict  the  melting  point  by  100  K. 


Fig,  13:  Density  versus  temperature 
for  NaC.1  at  P=1  bar.  The  dotted  line  is 
the  experimental  melting  point,  while 
the  simulation  results  (filled  circles) 
predict  a  melting  point  that  is  roughly 
100  K  too  high.  This  demonstrates  the 
problem  with  gradual  heating 
approaches  for  melting  point 
prediction. 


To  compute  the  true  thermodynamic  melting  point,  the  free  energy  of  the  liquid  and  solid 
phases  are  computed  with  respect  to  separate  reference  states  using  histogram 
reweighting  techniques5.  An  example  of  these  free  enrgy  curves  for  NaCI  is  shown  in  Fig. 
14.  Notice  that  each  curve  is  computed  with  respect  to  an  arbitrary  reference  state  in  its 
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own  phase  at  T-1080. 

Fig.  14:  Free  energy  curves  for  the  solid 
and  liquid  phases  computed  from 
histogram  reweighting  simulations 
performed  in  the  iso  thermal-  iso  baric 
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ensemble.  Note  that  each  curve  is  referenced  to  its  own  phase  at  1080  K. 

To  compute  the  melting  point,  the  two  curves  in  Fig.  !4  must  be  referenced  to  each  other. 
This  is  accomplished  using  a  three-step  thermodynamic  integration  path,  shown 
schematically  in  Fig.  15.  Once  the  free  energy  difference  between  the  solid  and  liquid 
curves  is  known  at  a  given  statepoint,  the  two  curves  in  Fig.  14  can  be  shifted  and  the 
intersection  will  yield  the  melting  point.  The  free  energy  difference  between  the  solid  and 
liquid  phase  for  NaCI  obtained  from  this  calculation  is  shown  in  Fig.  16. 


Fig.  15:  Three-stage  pathway  for 
computing  the  free  energy  between  the 
solid  and  liquid  phase.  In  (a),  the  liquid 
is  changed  to  a  weakly  interacting 
pseudo-supercritical  fluid.  In  (b)  this 
fluid  is  compressed  to  the  solid  density 
and  Gaussian  potential  wells  are  turned 
on  to  include  the  underlying  crystalline 
order.  In  (c),  the  Gaussian  wells  are 
gradually  turned  off  while  the 
intermolecular  potential  is  turned  back 
on,  yielding  an  ordered  crystalline  solid. 


Fig.  16:  Computed free  energy 
difference  between  crystal  and  liquid 
of  NaCI  at  P=1  bar.  Location  where 
AG  =  0  is  the  melting  point. 


The  P=1  bar  melting  point  had  been  previously  computed  for  this  model  by  Frenkel  et 
al.16.  Our  results  agree  quite  well  with  theirs.  The  high  pressure  melting  point  for  this 
model  had  not  been  computed  up  until  the  present  study. 


Given  a  single  melting  point,  the  Gibbs-Duhem  integration  procedure  can  be  used  to 
obtain  the  full  solid-liquid  coexistence  line,  This  was  done  for  NaCI  up  to  P=20.000  bar. 
Results  are  shown  in  Fig.  17  and  compared  with  experimental  data.  It  can  be  seen  that, 
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while  the  melting  point  at  P=l  bar  is  in  excellent  agreement  with  experiment,  the  mode! 
tends  to  overpredict  the  melting  point  at  higher  pressures.  This  is  not  surprising,  as  the 
model  was  parameterized  at  atmospheric  pressure.  It  points  out  both  the  power  of  the 
method  and  the  importance  of  considering  forcefield  parameterization  over  a  range  of 
conditions. 


Fig  7  7:  Full  coexistence  curve 
computed  for  NaCl  using  Gibbs - 
Duhem  integration,  along  with 
experimental  data.  It  can  be  seen 
that  the  model  tends  to  overpredict 
the  melting  point. 


The  method  was  extended  to  compute  the  melting  points  of  multi-atom  molecular 
systems.  The  molecules  1-H- 1,2, 4-triazole  and  benzene  were  examined  as  test  systems. 
The  former  molecule  is  a  precursor  of  energetic  ionic  liquids  and  crystal  structures  and 
melting  points  arc  known.  The  latter  has  a  highly  optimized  force  field  and  previous 
melting  point  simulations  using  non-free  energy  based  methods  have  shown  that  the 
computed  melting  point  does  not  agree  well  with  experiment. 


To  simulate  these  systems,  each  individual  atom  had  to  be  tethered  to  a  lattice  site.  This 


Figure  1 8:  Computed  probability  distribution 
of  the  triazole  center  of  mass  about  its  lattice 
site,  and  the  distribution  obtained  using  a 
harmonic  tethering  potential 


was  done  by  carrying  out  a  solid 
phase  simulation  and  then  Fitting 
harmonic  force  constants  for  each 
tether  that  reproduced  the  natural 
vibrational  motion  of  the  solid 
phase.  Fig.  18  shows  an  example  of 
the  actual  deviation  of  a  triazole 
molecule  from  its  lattice  site  during  a 
solid  simulation  and  the  fitted 
tethering  potential. 

The  tri azole  model  is  found  to  do  a 
poor  job  reproducing  the 

experimental  crystal  structure  and 
subsequently  yields  a  melting  point 
that  is  100  K  too  low.  The  benzene 
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model,  however,  was  chosen  based  on  ability  to  reproduce  experimental  data,  and  the 
resulting  melting  point  is  within  20  K.  of  the  experimental  value.  This  compares  well  with 
a  recent  study  that  found  a  much  larger  discrepancy  between  simulated  and  experimental 
melting  point  for  the  same  model  based  on  a  method  without  a  free  energy  analysis. 
Based  on  the  different  results  for  triazole  and  benzene,  the  sensitive  relationship  between 
force  field  and  crystal  properties  is  discussed  and  related  to  these  results. 

3.5.3  Relationship  to  the  Goals  of  the  Project 

This  work  has  resulted  in  a  significant  new  simulation  method  that  can  be  applied  to 
rigorously  compute  melting  points  of  molecular  systems.  Such  a  technique  is  directly 
applicable  to  the  ionic  liquids  research,  as  rapid  a  priori  prediction  of  melting  points  will 
greatly  assist  experimental  efforts  directed  at  discovering  low  melting  energetic  ionic 
liquids.  Because  the  method  is  rigorous,  it  allows  for  unambiguous  determination  of  the 
melting  point  and  hence  stringent  tests  of  intermolecular  force  fields.  We  are  currently 
using  the  method  to  simulate  melting  points  for  ionic  liquids. 
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